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Abstract 

We propose a new coarse-grained model for the description of liquid-vapor phase separation of 
colloid-polymer mixtures. The hard-sphere repulsion between colloids and between colloids and 
polymers, which is used in the well-known Asakura-Oosawa (AO) model, is replaced by Weeks- 
Chandler- Anderson potentials. Similarly, a soft potential of height comparable to thermal energy 
is used for the polymer-polymer interaction, rather than treating polymers as ideal gas particles. It 
is shown by grand-canonical Monte Carlo simulations that this model leads to a coexistence curve 
that almost coincides with that of the AO model and the Ising critical behavior of static quantities 
is reproduced. Then the main advantage of the model is exploited — its suitability for Molecular 
Dynamics simulations — to study the dynamics of mean square displacements of the particles, 
transport coefficients such as the self-diffusion and interdiffusion coefficients, and dynamic structure 
factors. While the self-diffusion of polymers increases slightly when the critical point is approached, 
the self-diffusion of colloids decreases and at criticality the colloid self-diffusion coefficient is about 
a factor of 10 smaller than that of the polymers. Critical slowing down of interdiffusion is observed, 
which is qualitatively similar to symmetric binary Lennard- Jones mixtures, for which no dynamic 
asymmetry of self-diffusion coefficients occurs. 

PACS numbers: 
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I. INTRODUCTION 



In the last few decades colloidal dispersions have been studied intensively as model sys- 
tems for the structure and phase behavior of fluids and solids. The large size of the colloidal 
particles allows for additional experimental techniques which are not applicable for atom- 
istic or molecular systems. Moreover, the colloid-colloid interactions can be "tuned" to a 
large extenti^^. For example, individual colloidal particles can be tracked through space 
in real time using confocal microscopy^. In colloid-polymer mixtures, where the depletion 
attraction between the colloids caused by the polymers 7 - 1 ^ can lead to a liquid-vapor type 
phase separatior*i&ii>i2, statics and dynamics of capillary wave-type interfacial fluctuations 
can be observed in real spaced. Wetting layers of the walls of containers can be stud- 
ied in detail^^£, and critical fluctuations can also be seen directly in optical microscope 
observations 1 ™. Very interesting nonequilibrium studies are also possible, such as shear- 
induced narrowing of interfacial widths^ and studies of spinodal decomposition 19 . 

In view of this wealth of experimental data on static and dynamic behavior relat- 
ing to liquid-vapor type phase separation in colloid-polymer mixtures, it is also de- 
sirable to provide a detailed theoretical understanding of these phenomena. In fact, 
many static aspects (including the understanding of the phase diagram and bulk critical 
behavior— >2i>22, interfacial fluctuations^ and interface localization transitions^ 1 ^, capillary 
condensat ion/e vapor at ion2^^'2Ii^i2£i^ and wetting^i^^^) can all be understood by the 
simple Asakura-Oosawa (AO) 7,8,9 model, at least qualitatively. In this model colloids and 
polymers are described as spheres of radius R c and R p , respectively. While there is a hard 
core interaction of the colloids both among each other and also with the polymers, the 
polymer-polymer interaction is assumed to be strictly zero. Thus, a suspension without any 
colloids but only polymers is just treated as an ideal gas of point particles which are located 
at the center of mass of the polymer coils. 

This model is very attractive due to its simplicity. It allows for various ele- 
gant analytical approximations^^^ 1 ^ 1 ^ as well as for efficient Monte Carlo simulation 
techniques^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^ 1 ^. However, the assumption that polymers do not inter- 
act with each other at all makes the AO model unsuitable for studying dynamical aspects 
of colloid-polymer mixtures. Thus, a different model is required to complement the corre- 
sponding very interesting experiments mentioned above^^^ 9 -. 
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Occasionally, computer simulations have been performed where the polymers were mod- 
eled explicitly as chain molecules either on the lattice^^ 1 ^ or as bead-spring-type chains 
in the continuum^. In general, these models are restricted to rather small chain lengths 
in order to keep the numerical effort manageable. In addition, only particle sizes in the 
nanometer range can be treated. However, one can use these simulations^ 7 - 1 ^ 8 - to justify an 
effective interaction between two polymer coils. Thus, polymers are described as soft par- 
ticles which can "sit on top of each other", but not without energy cost. The usefulness of 
such an effective potential has been amply demonstrated^I^^ML 1 ^. 

This consideration is the motivation for the present study. We define a model (Sec. [II]) 
which has a soft interaction potential between polymers, too, and is particularly convenient 
for both Monte Carlo^^ and Molecular Dynamics^ 1 ^ simulations. In Sec. II I II the static 
properties of the model are evaluated and compared to corresponding results^ 1 ^ for the 
standard AO model 7,8,9 . Section [TV] presents our data for the mean square displacements 
of the particles as well as intermediate scattering functions. We also discuss the resulting 
self-diffusion and interdiffusion coefficients while Sec. [V] summarizes our conclusions. 

II. A SOFT VARIANT OF THE AO MODEL 

A potential of type U(r) = Uq exp[— (r/R g ) 2 ] describes the effective interaction between 
two polymer coils in dilute solution under good solvent conditions. This result can be 
obtained by calculating the partition function of the two chains under the constraint that 
the distance r between the centers of mass of the coils is fixed. The prefactor Uq is of the 
order of the thermal energy^^^^i and R g is the radius of gyration of the chains. Similarly, 
the interaction potential between a polymer chain and a colloidal particle can be obtained. 

However, the situation becomes slightly more involved at higher polymer concentrations 
where many coils overlap and the temperature of the polymer solution (in comparison with 
the Theta temperature^) also plays a role. Then it is no longer possible to give a simple 
explicit description for the polymer-polymer interaction from first principles. Additionally, 
it is more convenient for computer simulations to have a potential which is strictly zero if r 
exceeds some cutoff r c . Therfore, we did not use any of the approximated effective potentials 
derived in the analytical work— li&i&il Instead we chose a potential that has qualitatively 
similar properties, but is optimal for our simulation purposes. For the colloid-colloid and 
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colloid-polymer potential we took the Weeks- Chandler- Anderson (WCA) potential 47 , mod- 
ified by a smoothing function S 



with 
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Here, e a p controls the strength and a a p the range of the (repulsive) interaction potential 
which becomes zero at r c>a p and stays identically zero for r > r CjQ/ g with r CjQ/ 3 = 2 1//6 er Qj g. Fol- 
lowing previous work on the AO model^^^ 1 ^ 1 ^ 4 - 1 ^ 1 ^ 1 ^, we chose a size ratio q = cr pp /cr cc = 
0.8 between polymers and colloidal particles and 



(T C p = 0.5(cx cc + (T pp ) = 0.9cr cc . 



(3) 



The parameter h of the smoothing function is taken as h = 10 cr cc and e 



- cp 



1. In 



the following, we choose units such that k^T = 1 and a cc = 1. Note that the smoothing 
function is needed in Eq. (CQ) such that U a p{r) becomes twofold differentiable at r C)CC and 
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without affecting the potential significantly for distances that are not very close to 



the cutoffs. Without S the force would not be differentiable at the cutoff distances and 
hence a noticeable violation of energy conservation would result in microcanonical Molecular 
Dynamics (MD) runs^ 4 -^. 

For the soft polymer-polymer potential the following somewhat arbitrary but convenient 
choices are made: 
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Note that the expansion in the square bracket of Eq. (jlj) is essentially a polynomial fit to a 
cosine function, which is shifted by unity and the angle of which varies from to 7r when r 
increases from zero to ^ c ,pp- However, while the cosine function is smoothly differentiable only 



once at r = and r = r CjPP , Eq. (J4]) is twofold different iable. Of course, U pp {r > r CjPP ) = 0. 
Note that the choice (jSJ) differs from the original AO model only by replacing the original 
hard core interactions U cc (r) and U cp (r) by smooth interactions, Eqs. ([I])-©, while polymers 
are still strictly non- interacting. For the choice (0), which is the choice used for the MD 
work, the energy varies from U pp (r = 0) = 1/2 k^T to zero, Fig. [TJ With these choices of 
potentials the application of MD is straightforward and efficient, but also the application of 
grand-canonical Monte Carlo methods is still well feasible. 

However, as in our earlier study of static and dynamic critical phenomena of a symmetri- 
cal binary Lennard- Jones mixture^^^L it is advantageous to apply Monte Carlo methods 
(in the fully grand-canonical (iV/i c /z p T)-ensemble, with p c , p p being the chemical potentials 
of colloids and polymers, respectively, in the present case) to determine the static phase 
diagram of the colloid-polymer mixture. In particular, determining the critical densities of 
colloids p£ nt and polymers p£ nt (where densities are defined in terms of the particle numbers 
of colloids N c and polymers N p in the standard way p c = N c /V, p p = N p /V, V being the 
volume of the simulation box) is a nontrivial matter. In the context of MD simulation, start- 
ing systems at states which fall within the two-phase coexistence region and beginning with 
an initially homogeneous distribution of both types of particles lead to a phase separation 
into a "liquid-like" phase (with densities p c , pp) and a "vapor-like" phase (with densities p", 
Pp). Of course, a priori all the values of densities along the coexistence curve in the (p c , p p )- 
plane are unknown and simulating the dynamics of spinodal decomposition is a complicated 
and notoriously slow process^ 1 ^. Moreover, when approaching the critical point (from the 
one-phase region or along the coexistence curve) simulations in the canonical ensemble suffer 
severely from critical slowing down^ as discussed 

Thus, it is very desirable to study the phase behavior by Monte Carlo simulations 
in the grand-canonical ensemble, which turned out to be very useful for both the sym- 
metrical binary Lennard- Jones mixture^ and the standard AO model^ 1 ^. a number 
of well-equilibrated system configurations as initial states for strictly microcanonical MD 
runs^ 1 ^ 1 ^ 1 ^, one realizes averages corresponding to a well-defined temperature T without 
the need to augment the MD code by a thermostats^. However, already for the standard 
AO model straightforward particle insertion Monte Carlo moves, which are necessary to 
realize the grand-canonical ensemble, are almost always rejected due to the large density of 
polymers p p in the system^ 1 ^. For the standard AO model Vink and Horbach^ 1 ^ could 
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cope with this difficulty by implementing a cluster move. A similar cluster move is used 
here (Fig. [2]) . To maximize the efficiency of this algorithm, always all polymer particles are 
removed in the depletion zone when a colloid is inserted. The radius of the depletion zone 
was cr cc + r C)Cp . At most m — 10 polymers would be inserted or removed in one attempted 
cluster move. 

When this algorithm is applied to soft potentials, slight modifications of the implementa- 
tion are required: Colloid deletion attempts must always be rejected if any center of polymer 
particles is located in the depletion zone. Otherwise, colloid insertion and deletion moves 
are no longer symmetric and detailed balance is violated. Note that this problem does not 
occur in the original AO model because polymer particles can never "overlap" with colloid 
particles. Note that the algorithm is still ergodic, because "overlaps" with the colloidal 
particles can be obtained by removing adjacent colloids and filling the void with polymer 
particles. Nevertheless it is still recommendable to mix cluster moves with local moves. 

We wish to compare our results with the original AO model (with hard core interactions), 
where it is standard practice to use the packing fractions r) c , r] p as variables, 

Vc = pcV c , r] p = p p V p , (7) 

where V c and V p are the volumes occupied by a colloid and polymer, respectively, with 
V c = 7tg^ c /6 and V p = Trd pp /6, d cc and d pp being the diameters of colloids and polymers. For 
this comparison it is hence useful to define an effective diameter of colloids and polymers of 
our model using the approach of Barker and Henderson^ 
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Using Eqs. ©, © in Eq. (JED yields 

d cc = 1.01557cr cc (9) 

and d cp = 0.9d cc . We also use d pp = 0.8d cc and consequently derive the following formulas 
to convert our densities into packing fractions 

r] c = 0.54844cx c 3 c p c , ?7 P = 0.28080(J c 3 c p p . (10) 

Hence, the polymer reservoir packing fractior*iii2Si>2i j s gi ven by ^ — V p exp(fi p /kBT) = 
0.28080 al c exp(/i p /fc B T) in terms of the chemical potential fi p of the polymers. Of course, 
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for interacting polymers the notion of rf^ loses its original meaning, but we continue to use 
?7p as defined here for the sake of comparability with the standard AO model. 

The technical aspects of grand-canonical Monte-Carlo simulations of phase equilibria 
and critical phenomena in colloid-polymer mixtures have been described in detail in the 
literature^^i^ 1 ^. Therefore, we recall only very briefly the most salient features. Choosing 
/i p and hence 77^ as a parameter, the chemical potential /i c of the colloids is varied and the 
distribution P(i] c ) of the colloid volume fraction is sampled, applying the cluster algorithm 
mentioned above. For rf p sufficiently less than the critical volume rf p cvlt (note that rj^ plays 
the role of inverse temperature, when the phase diagram of the colloid-polymer mixture is 
compared to the vapor-liquid phase separation of a molecular system) P(i] c ) has a single 
peak at (r) c ) and the task is straightforward. For 77I > r£ crit and \i near fi coex , however, P(rj c ) 
is a doubly-peaked function where (apart from finite size effects^^i^ 1 ^ 1 ^ 1 ^) the positions of 
the two peaks correspond to the volume fractions of the vapor-like phases at the coexistence 
curve, r^{vQ and ff c {rQ. Since the distribution P(r] c ) develops a very deep minimum in 
between these peaks^ 1 ^, sampling is, however, not completely straightforward. An efficient 
way to overcome this difficulty is provided by successive umbrella sampling 60 , which was 
applied in this work. The chemical potential /x coex at vapor-liquid coexistence is then given 
by the equal weight ruleQ, i.e. the areas underneath the peaks corresponding to the vapor- 
like phase and the liquid-like phase have to be equal. The order parameter m of the phase 
transition can then be identified as 

m = 7] c - (r} c ) , (11) 

with (i] c ) being the average of P{f] c ) including both peaks at /U coex - A convenient tool to find 
the critical value rj^ crit is based on the analysis of moment ratios M, U defined as^>^ 

M= « U=M (12) 
{\m\) 2 {m z ) z 

while following a path along fi coe x{f] r p ) for different linear dimensions L of the cubic simulation 
box. The critical value ^p Crit is determined found from the intersection of these curves. 
Figure [3] gives an example for the present model. Note that the procedure described above 
is also operationally well defined for a range of values vL < ^p Crit , since due to finite size 
effects, the distribution P(r] c ) is double-peaked over some range in the one-phase region as 
we }j56i58 jj. can j-jg recognized from Fig. [3] that a rather well-defined intersection point occurs 
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for r)p Crit = 1.282±0.002. However, this intersection does not occur at the theoretical valued 
M ~ 1.239 but at a somewhat lower value M c g « 1.21. This discrepancy is due to various 
corrections to finite size scaling, in particular the so-called field mixing effects^^. In the case 
of the standard AO model, a very similar discrepancy occurs as wei l 20 ' 21 . Since in the latter 
model no potential energy is present, the field mixing does not involve a coupling between 
energy density and density as for ordinary fluids 63,64 but rather a coupling between colloid 
density and polymer density. So the order parameter (in the sense of a scaling fieldM-^) is 
in a strict sense not given by r] c alone (as assumed in Eq. (fTT!) ). Instead, a suitable linear 
combination of r] c and r] p needs to be constructed. However, we have not done this in the 
context of finding the critical point since there is ample evidence in various systems 6 ^* 6 ^ 6 ^ 
that the simple cumulant intersection method as illustrated in Fig. [3] does yield the critical 
point with a relative accuracy of a few parts in a thousand, which suffices for the present 
purposes. 

III. STATIC PROPERTIES OF THE SOFT VERSION OF THE AO MODEL 

As discussed in Sec. [TTJ the first step of the Monte Carlo study consists of the estimation 
of the coexistence curve and the critical point. For the two models defined in Eqs. ((SD,® 
we found for e pp = 

<crit = 0.760, Vc , crit = 0.136, 77 P)Crit = 0.354 (13) 

and for e pp = 0.0625 

7£ iCrit = 1.282, ry C)Crit = 0.150, r) P;Crit = 0.328 . (14) 

Since the accuracy of these numbers is about ±0.002, we conclude that model (i), the "soft 
AO model" , is within our errors not distinguishable from the original AO model with hard 
core interactions for which the analogous results are 2 ^ 2 ^- 

V;, crit = 0-766, 77 CiCrit = 0.134, ry p , crit = 0.356 . (15) 

This coincidence between the soft AO model and its hard core version is also seen in the 
coexistence curve, which is compared in reservoir representation in Fig. HI The coexistence 
curve of the model with interacting polymers is substantially different, of course, as expected 
from Eq. (114]) . 
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However, when we study phase coexistence as a function of all experimentally accessible 
variables f] p ,f] c , differences between the three models are rather minor (Fig. [5]). It appears 
that near criticality the main effect of "switching on" the polymer-polymer interaction is to 
shift the critical point along the coexistence curve of the AO model to the higher value of 
^c,crit (and correspondingly lower value of 77 PiCr it mentioned in Eq. (j!4p . Further away from 
the critical point the coexistence curve of the interacting polymer model predicts somewhat 
lower polymer packing fractions along the "vapor branch" and somewhat higher polymer 
packing fractions along the "liquid branch". The result that the critical packing fraction 
of polymers is about twice that of the colloids is similar to what was observed in a recent 
experiment^. Note however, that in these experiments a size ratio of polymers to colloids 
of q = 1.04 (rather than q = 0.8) was employed which affects t/ c cv h 

(?7c,crit ~ 0.10 was found 

inlft). 

In Fig. [5] we have also indicated the state points at which micro-canonical MD runs 
took place. We have fixed the number of colloids N c in a volume (of size 27 3 ) such that 
Vc = ??c,crit ~ 0.15 (which corresponds to iV c = 5373). Then the number of polymers was 
varied from zero up to iVp = 22734. The MD runs were carried out with the Velocity Verlet 
algorithm^ 1 ^ and a time step 5t = 0.0005(cr^ c m c /e cc ) 1 / 2 . The masses of colloids and poly- 
mers are equal and units of time are chosen such that m c = m p = 1. In the production 
runs used for the computation of time-displaced correlation functions no thermostat was 
applied, so a microcanonical ensemble respecting all conservation laws applies. Starting 
configurations were generated as follows: First, a random configuration was generated in 
a box of linear dimension L = 9 and periodic boundary conditions. The system was equi- 
librated at T = 1 for 20 million time steps with a simple velocity rescaling according to 
the Maxwell-Boltzmann distribution. Then, the system is enlarged from L = 9 to L = 27 
by replicating it three times in all spatial directions. Now periodic boundary condition for 
L = 27 only are applied. Equilibration is continued for 2 million time steps, again with 
a Maxwell-Boltzmann thermostat. During this equilibration, the original periodicity with 
L = 9 is quickly lost. The production runs for static averages are done without applying any 
thermostat. First, 5 million time steps are performed during which (at eight different times) 
statistically independent configurations are stored. These serve as starting configurations 
for eight independent simulation runs, each with 5 million steps, for the computation of 
static averages. During each run, 500 configurations are analyzed in regular intervals. Thus, 
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4000 statistically independent configurations are averaged over for the computation of the 
structure factor. 

From now on, we denote colloids as A-particles and polymers as B-particles. For all 
simulated state points the partial structure factors were computed, 

S a p = ( ^^exp(zg- fij) ) , aeA,B, (16) 
\i=i j=i I 

with N = Na+Nb- These results, presented in Figs. [6(a)- (c), show that the partial structure 
factor for colloids (Fig. [6(a)) displays an oscillatory structure with a first peak near q w 6.5, 
which is typical for the packing of hard particles in a moderately dense liquid. The partial 
structure factor due to polymers (Fig. [6(c)) exhibits much less structure in the range of 
large q as expected, since for the potential, Eq. (j4j), the polymers can still overlap rather 
easily. All these partial structure factors show a strong enhancement at small q, reflecting 
the critical scattering due to the unmixing tendency between colloids and polymers when 
the critical point is approached. Note that the partial structure factor due to interference 
of the scattering from colloids and polymers (Fig. [6](b)) also shows oscillations at large q, as 
does the scattering from colloids alone (Fig. E[a)). 

From the partial structure factors it is useful to construct combinations that single out 
number-density fluctuations Snn(i) and concentration fluctuations Scciq), defined via^ 
(x A = N A /[N A + N B },x B = l- x A ) 

Snn(q) = Saa(?) + 2<Sab(<?) + S BB (q) , (17) 
Scc(q) = xlS AA (q) + x 2 A S BB (q) - 2x A x B S AB (q) . (18) 

In addition, it is of interest to consider a structure factor relating to the coherent interference 
of number density and concentration fluctuations 68 , 

<SVc(g) = x B S AA (q) - x A S BB (q) + (x B - x A )S AB (q) . (19) 

Figure [7] shows that all three structure factors show a strong increase at small q, reflecting 
the critical scattering as the critical point is approached. Additionally, at large q they 
display oscillations. The behavior seen in Fig. [7| differs very much from the behavior found 
for the unmixing of the symmetric binary Lennard- Jones mixture^^^i. In the latter 
case Snn(q) was not sensitive to the critical fluctuations at all, which showed up in Scc{q) 



11 



only. Likewise, Scciq) was insensitive to the way how the particles are "packed" in the 
liquid, i.e. there was no structure at large q. In addition, almost no interference between the 
scattering from concentration and density fluctuations could be seen. Hence, Snc(q) was 
very small, while in the present model Sccil) an d Snc{q) are of the same order of magnitude. 
These observations clearly show that neither the total density in the system, nor the relative 
concentration of one species is a "good" order parameter of the phase separation that occurs. 
(Likewise, Fig. [6] shows that neither the colloid density alone nor the polymer density alone 
are "good" order parameters since both densities reflect the critical scaling in a similar way.) 
Of course, from the phase diagram (Fig. [5]) such a problem is expected since the shape of 
the coexistence curve shows that the order parameter is a nontrivial linear combination of 
both particle numbers N p , N c . 

In order to deal with this problem, we introduce a symmetrical matrix formed from the 
structure factors Saa(<?), Sab(q) an d S BB (q) 

m = ( S «<») «»<») | (20 ) 
\Sab(?) S BB (q) 

and diagonalize this matrix to obtain its diagonal form 

5«M=f S+(?) ° ], (21) 

with 



S±(q) = \[S A A(q) + S BB (q)] ± ^\[S AA (q) - S BB (g)] 2 + S 2 AB (q) . (22) 

Figure M shows a plot of S+(q) and <SL(?) versus q. This plot shows that this procedure 
indeed resulted in a decoupling of the order parameter fluctuations (which show a critical 
enhancement as q — > 0), as being measured by S + (q), and the noncritical "particle packing" 
fluctuations, measured by S-(q), which show the characteristic oscillatory structure of a 
noncritical fluid. In the case of the symmetrical LJ mixture the transformation from the 
number density fluctuations of A and B particles to the structure factors measuring the 
fluctuations of the total density of particles and of their relative concentrations is unam- 
biguous. In the case of the colloid-polymer mixture it is none of these variables which plays 
the role of an order parameter, but a different linear combination of both local densities of 
A and B particles, related to the eigenvector corresponding to S + (q). We can give this fact 
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a plausible interpretation by constructing two linear combinations of the operators p\{q), 
Pb(q), defined via p a (o) = Si=i ex P( i( f' as follows 

^j(q) = ap A (q) +bp B (q) , (23) 

<P(q) = a'pA(q)+b'p B (q). (24) 

The coefficients a, b are defined such that at the critical point the densities lie tangential 
to the coexistence curves. Coefficients a', b' are chosen such that the densities vary in a 
perpendicular direction to this slope. When we construct the structure factors (Fig. [9]) 

SW(?) = ^mq)\ 2 ) , S„{q) = ^<|0(g)| 2 > , (25) 

one recognizes that S^(q) is very similar to S + (q) and S^q) very similar to S-(q). The 
structure factors defined in this manner are not strictly identical to S+(q), S-(q). With 
increasing distance from criticality the relative weights b/a, b'/a' of the components of the 
"order parameter components" ip{q),4>{q) change. 

For q — > all those structure factors that show a critical increase can be described by 
the well-known Ornstein-Zernike behavior. This is illustrated in Fig. [TOj as an example, for 
the concentration, fitting 1/Scc{q) versus q 2 at small enough q (q 2 <C 2) to the relation^i 

Sccti) = {k B T X ccY l [l + q 2 e cc + •••]> 9^0. (26) 

Here, \cc is the "susceptibility" describing the magnitude of concentration fluctuations and 
^cc the correlation length. The various susceptibilities relating to the various structure 
factors defined above and the associated correlation ranges are shown in Fig. [TTJ It is 
gratifying to note that indeed the "susceptibility" related to S+(q) is the largest susceptibility 
that can be found, while the estimates for the correlation lengths are all equal (within 
statistical errors). Due to the coupling between variables, there is only a single correlation 
length in the problem. 

In Fig. [11] we have included two theoretical predictions in the log-log plot for the crit- 
ical exponents, one is a slope corresponding to the standard Ising exponents (that are 
observed in the grand-canonical ensemble, where only intensive thermodynamic variables 
are held constant). The other slope shows the exponents if "Fisher renormalization" 
occurs. To remind the reader of this phenomenon we note that the response function 
X = (dN c /dp)T, f i p \[i=[ JjcIit /N that is observed via Monte Carlo from the fluctuation relation 

Wx = N~\(N 2 ) - (N c ) 2 ) Tt , p (27) 
13 



differs from xcc since fluctuations differ in different ensembles of statistical mechanics. 
While xcc was estimated from Eq. fl26|) which refers to the ensemble where N p = const, 
Eq. (I2"7j) refers to the ensemble where /i p = const. Sinco^ 

N p = iVp >crit + a(/i p - yUp^it) 1- " + b(/i p - /i p ,crit) + • • • , (28) 

where a is the specific heat exponent and a, b are constants. Very close to the critical point 
we have a singular relation between (N p ) — iVp )Crit and fi p — ji p iCrit , namely 



N p / /i p 



l-a 



1 oc - 1 . (29) 

i*p,erit V/ipjCrit / 

Therefore the power laws in the grand-canonical ensemble^^i 

X oc f-^ - 1 V , £ oc f-^ - l) " (30) 

translate into power laws with "Fisher renormalized"— exponents in the microcanonical 
ensembles where N p , N c are constant 

Xcc oc e-^ 1 -") , i cc oc e^ 1 ^ . (31) 

However, since the regular third term on the right hand side of Eq. (|28|) is comparable to the 
(singular) second term that was only used in Eq. (129]) . except if one works extremely close 
to yU Pi crit, it is difficult to ascertain whether or not the simulation data shows any signature 
of Fisher renormalization. High precision simulations for very much larger systems would 
be required to clearly resolve this issue. This task, however, is not possible with presently 
available computer resources. 

It is also useful to recall that susceptibilities observed in the grandcanonical ensemble 
differ from those extracted from structure factors in the canonical ensemble. To interpret 
this difference, we start from the standard relation for the grandcanonical partition function 

Z gc (p, /ip, V, T) = J2 (f^) £ exp (fffil Z C (N C , N p , V, T) (32) 

N c =0 V B / iVp=0 V B / 

from which one straightforwardly derives the following fluctuation relations (p c = (N c )/V): 
kB T d J^L = (A?) _ W , f d J?A = {d) _ (pc)2 (33) 

fc ^ = W-W, (34) 
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and 



k V T ~^ 

dfip 
or 

k B Td( Pc ) 



k B T^^ = (N c N p )-(N c )(N p ) 



(Pcpp) - (Pc>(Pp) 



(35) 



V dfj,p 

It is this mixed susceptibility describing the correlations between the fluctuations of colloid 
and polymer number which enters the difference between the susceptibilities in the canonical 
and grandcanonical ensemble. A simple calculation yields 



con _ y ( 9(pc) 

at, Pp — Ar 



N V dfi 



X c t oU 



V 
N 



d( Pc ) 
. d Pp . 



' d( Pp ) 
dfi p 



(36) 



Similarly, 



poi = y_ ( d(p P ) 

V 
N 



X.T,n 



T,pc 

d(p c ) 
dp P 



' d(Pc) 



(37) 



In fully grand-canonical simulations, as we have carried out in the present work, it is 
possible to extract all susceptibilities of interest from a study of the joint distribution function 
P(N C , N p ). In the one phase region and for large enough linear dimensions L this function 
is a bivariate Gaussian in the variables N c — (N c ), N p — (N p ), see Fig. [12] for an explicit 
example. Then the fluctuations Xtp p anc ^ Xt°p c can ^ e extracted from the half-widths of 
these distributions along the abscissa direction (p p = N p /V = const) and ordinate direction 
(p c = N c /V = const), respectively. For the grand-canonical simulations described in Fig. 12, 
for example, we obtain XTp p ={p p } = 0-047 and 0.05 for the second equation in (36) thus 



confirming our calculations. Figure [12] illustrates again that none of these susceptibilities 
should be regarded as the order parameter susceptibility x+ : rather the latter is the half- 
width along the main axis of the ellipsoidal contours P{N C , N p ) = const in Fig. [12] 
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IV. DYNAMICS OF COLLOID-POLYMER MIXTURES 



From the MD runs it is straightforward to obtain the incoherent intermediate scattering 
functions F"(q,t) defined as (a = A, B) 

KM = ^ 5>xpHg • m) - f 4 (0)])> (38) 
as well as time-displaced mean square displacements of the particles 

= at £<fu*) - r ~^(°)] 2 > • ( 39 ) 

In the MD framework the average (• • • ) stands for an average over the origins of time, t = 
(we have used 8 statistically independent runs and two time origins per run, thus we have 
averaged over 16 time origins). Figure [TBI shows typical data for both small and large q. 
For the colloids there is some uniform slowing down of F°(q,t) at small q as N p increases, 
while for large q (near the first peak of S a p(q)) the decay occurs in two parts: the first 
part (for Fg(q,t) > 0.8) is basically independent of N p , while for F£(q,t) < 0.5 the curves 
distinctly splay out. In contrast, the analogous function for the polymers Ff(q,t) seems to 
be practically independent of N p , irrespective of q. 

A similar asymmetry between the dynamics of colloids and polymers is also seen in the 
mean square displacements. Since we expect for large times the Einstein relation to hold, 

g a (t) = QD a t, t^oo, (40) 

we analyze the derivative (1 / 6)dg a (t) / dt (Fig. [TJJ. From the plateau of this quantity at large 
times, one can see that g a (t) approaches its asymptotic behavior for colloids monotonically 
while for polymers there is an overshoot for intermediate times, 1 < t < 10. In the regime 
of this transient maximum the data depends rather distinctly on N p . In the asymptotic 
regime (t — > oo) the dependence is much weaker. The time range where this overshoot 
occurs is related to the crossover from ballistic to diffusive motion. For i < 1 both colloids 
and polymers show a ballistic behavior, g a oc t 2 , as expected^ 1 ^. Of course, no such 
behavior is expected for real colloid-polymer mixtures, where the solvent molecules (no 
explicit solvent is included in our simulations, of course) damp out the "free flight" motion 
present in our model. Instead, one would find another diffusive motion controlled by the 
solvent viscosity. Figure [15] shows that the resulting selfdiffusion constants are of similar 
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magnitude for small N p (very far from iV P)Cr i t ) but differ by almost an order of magnitude 
when N PtCrit is approached. 

Finally, we consider the interdiffusion between colloids and polymers. Defining the center 
of mass coordinate of the particles of species a as R a (t), we note that interdiffusion is related 
to the following mean square displacement (m A = itib)^ 

^int(t) = <[rint(t)-r int (0)] 2 ) 
N A \ 2 N A N B 



1+ t N^^-^n- (41) 



Note that R A (t)-R A (0) is computed via the integral J * V A (t')dt' with V A {t) = N A X ^]f A v^t) 
the center of mass velocity of component A and Vi(t) the velocity of particle % at time t. In this 
manner one obtains the difference R A (t) — R A (0) in an origin independent representation^ 1 ^. 
The Onsager coefficient A relating to interdiffusion can be expressed as 

A = lim A(t) , A(t) = (6t)-^ int (t) . (42) 

t — >oo 

The interdiffusion constant -Dab, which describes how concentration fluctuations in the 
binary (A,B)-system relax, is then given as the ratio of the Onsager coefficient A and the 
"concentration susceptibility" Xcc, where 

_ x A (l - x A ) . _ S cc (q = 0) u „s 

VAB - j rj, A , XCC - 7—= . (4:6) 

KbTxcc k B T 
Note that theory^^^^^ predicts that A contains two terms, a background term A& which 
is nonsingular and stays finite at the critical point and a critical term AA which diverges at 
the critical point, 

'/p 



A = A 6 + AA , AA oc 1 ^- (44) 

V ^p,crit / 

with an exponent v\ ~ 0.567— iZ2j§£. In fact, a recent MD study of the critical dynamics 
of the symmetric binary Lennard- Jones mixture 5 ^' 51 yielded results compatible with this 
theoretical prediction, Eq. (j44j) . This allows to estimate the noncritical background term 
Afy at the critical point, too. Thus, it is also of great interest to study the behavior of A 
when we approach the critical point of our model (Fig. HT)) . Here, we have also included the 
simple prediction of the Darken equation^, 

A = x A D B + (1 - x A )D A . (45) 
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While very far from criticality (1 — ^ p /%,crit) > 0.6, Eq. (145]) indeed describes the simulation 
results accurately, it underestimates A strongly for r] p closer to % jCr it, and clearly Eq. (14"5]) 
violates Eq. (I44p . Thus, Darken's equation^ fails near the critical point of a fluid binary 
mixture as it was already noted for the binary Lennard- Jones mixture^. 

Thus, we see from Fig. [T7] that for our asymmetric mixture we also find evidence for 
a singular behavior of the Onsager coefficient for interdiffusion. However, the statistical 
accuracy of the data for A does not warrant an attempt to estimate the dynamic critical 
exponent v\ (in particular since this is rather difficult here to estimate A;,). The statistical 
effort invested is just enough to allow an approach of criticality up to about e = 1 — 
N p /N P!CTit ps 0.03, but not closer. In order to allow meaningful estimates of £, Xcc, an d A, 
the time r run of a simulation run must be at least about an order of magnitude longer than 
the time r needed for a concentration fluctuation to relax via interdiffusion. This time is 

r = (6D AB )~r = k ~^f- ■ (46) 

From Figs. fTTT a).(b) and [T71 we see for e = 0.03 that k-^Tx ~ 40, £ ps 6, and A ps 1, yielding 
r ps 240. Since r mn = 2500, the run at e = 0.03 is just long enough, but data closer to 
criticality cannot be used. The estimate Eq. (j46p is compatible with a direct examination 
of A(t), Fig. dni where we see that far away from criticality a plateau is only reached when 
r ps 100. 

Another condition for the validity of our result is that the initial periodicity width L- m i t = 9 
has fully relaxed. This equilibration time of our system is estimated in analogy to Eq. (146]) 
as T cq = (6-DAB) -1 -^nit ~ 540 for e = 0.03. The actual equilibration time of 10 3 MD time 
units indeed exceeds this estimate by a factor of about two. So our data should be valid but 
it is hardly possible to approach criticality closer. Finally, since no attempt of a finite size 
scaling analysis of the dynamical properties is made here (unlike^ 1 ^) , we have to require 
that L 3> 2£ at the states of interest. Though this condition holds for e = 0.03, it would fail 
if we approach the critical point much closer. From this discussion we see that a substantially 
larger computational effort would be required for a more detailed analysis of the dynamic 
critical behavior of this model. 
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V. SUMMARY 



In this paper, a model of colloid-polymer mixtures has been introduced and studied, 
which uses continuous potentials between all types of particles. Nevertheless, it still resem- 
bles closely the Asakura-Oosawa (AO) model, as far as static properties are concerned. The 
chosen potentials (Eqs. ([I])-©) are clearly somewhat arbitrary: the choice of these potentials 
was motivated by the desire that the model should be suitable for grand-canonical Monte 
Carlo methods to accurately establish the phase diagram (Figs. IU [5]). In addition, it should 
allow for a convenient and physically meaningful application of Molecular Dynamics tech- 
niques. In this way, both static and dynamic behavior of such a phase-separating strongly 
asymmetric binary mixture, where phase separation is mainly driven by entropic depletion 
effects, has become accessible to a computer simulation study. 

Previous simulation studies have mostly been concerned with the phase diagram of the 
AO model and related models as well as the interfacial tension between coexisting polymer- 
rich and colloid-rich phases. The present work contains a detailed analysis of the various 
static structure factors Saa{q), Sbb(<i) and Sab(<z) that one can define in such a binary (AB) 
mixture, and suitable linear combinations of them that single out the order parameter of the 
unmixing transition. Unlike other models of (almost incompressible) binary mixtures, in the 
present system the relative concentration of one species is not the proper order parameter. 
Instead the order parameter is a nontrivial combination of polymer density and colloid 
density fluctuations, which can be found from diagonalizing the structure factor matrix. 
From this analysis, we can study the onset of the critical divergence of both the order 
parameter "susceptibility" and correlation length. Roughly, these results are compatible 
with the expected Ising-like criticality. Fine details such as whether Fisher renormalization 
of critical exponents occur can, unfortunately, not be clarified, since our data are restricted 
to relative distances from the critical point exceeding 0.04. Of course, for reliable statements 
on critical exponents data somewhat closer to the critical point are indispensable, but at 
present not yet available. 

A central part of our study concerns the analysis of time-dependent quantities, interme- 
diate incoherent structure factors and mean square displacements, and their analysis. While 
the self-diffusion constant of the colloids is decreasing monotonously with increasing polymer 
density, surprisingly the self-diffusion constant of the polymers shows a slight increase. Thus, 
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there is a pronounced dynamic asymmetry of our model. The Onsager coefficient relating 
to interdiffusion is also obtained, and qualitative evidence for a critical divergence is found, 
thus invalidating the simple Darken equation for this system. However, the present data do 
not yet allow an accurate estimation of dynamic critical exponents for the colloid-polymer 
mixture. More efficient algorithms (or significantly faster computers) will be needed for a 
more definite study of critical behavior in our model system. Nevertheless, we hope that our 
study will motivate related experimental work on the dynamics of colloid-polymer mixtures, 
to which some of our findings could be compared directly. 
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FIG. 2: Schematic illustration of the cluster move applied for the grand-canonical Monte Carlo 
simulation. One colloidal particle (black sphere) is replaced by several overlapping polymer particles 
(dark grey) with randomly chosen center of mass positions inside the depletion zone of the colloid 
(indicated by the dashed circle). Other polymers in the environment are shown as light grey 
spheres. The double arrow indicates that the inverse move is implemented as well. 
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FIG. 3: Cumulant ratios M and U as a function of polymer reservoir packing fraction rjL for the 
model with interacting polymers (with 8e pp = 0.5). Three linear dimensions L are included (L 
is measured in units of cr cc ). The horizontal lines indicate the universal values^ M ~ 1.239 and 
U ~ 1.589 which M and U should acquire at criticality for every system in the universality class 
of the three-dimensional Ising model. 
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FIG. 4: Coexistence curves of the hard core AO model (fro m 20 ' 21 , full curve), the soft AO model 
(0), and the model with interacting polymers (0) in the plane of variable rj c , rj^ (reservoir repre- 
sentation). The open circle marks the locus of the critical point for the hard core and the soft AO 
model with e = 0. The full dot shows the critical point for the model with interacting polymers 
(e = 0.0625). 
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FIG. 5: Phase diagrams of colloid-polymer mixture models in the plane of variables r/ p (polymer 
packing fraction) and r/ c (colloid packing fraction). The original AO model (full curve), the soft AO 
model (standing crosses), and the model with interacting polymers (dashed curve) are compared. 
The triangles indicate state points at which MD runs were performed. Each coexistence simulation 
took 24 hours on a 32-core Power4 cluster (1.7 GHz). 
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FIG. 6: (a) Partial structure factor Saa(q) describing the scattering from colloids only, choosing 
ijc = f?c,crit = 0.150 and various choices for r] p as indicated for a simulation box of linear dimension 
L = 27. (b) Partial structure factor Sab(q) describing the interference in the scattering from 
colloids and polymers, (c) Partial structure factor Sbb(q) describing the scattering from polymers 
only. The values for the polymer packing fractions are rj p = 0.065, 0.129, 0.197, 0.255, 0.318 and are 
also used in Figs. ITlfTUl 
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FIG. 7: (a) Number density structure factor Snn{q) calculated from the data of of Fig. [H (b) 
Concentration structure factor Scc(q)- ( c ) Density-concentration interference structure factor 
Snc{q)- 
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FIG. 9: Structure factors S^(q), (a), and S^^q), (b), plotted versus q for the same values of r] p 
as in Fig. EJ Coefficients are a = -0.24, 6 = 0.97, a' = -0.97, b' = -0.24. 
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FIG. 11: Log-log plots of (a) k^Tx and (b) £ versus e = 1 — 7/p/r? p , C rit from MD simulations. Dashed 
and dashed-dotted lines indicate power law fits with (a) the exponents k-gTx oc e -7 and e _7r and 
(b) £ oc eT v or eT Vr , where 7 = 1.24 and v = 0.63 are the standard Ising exponent a 70 ' 71 while 
7 r = 7/(1 — a) and v T = u/(l — a) are the Fisher renormalized exponents— where a ~ 0.11 is the 
critical exponent of the specific heat^. 
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FIG. 12: Contour plot of a two dimensional probability distribution P(N C ,N P ) in the one phase 
region (e = 0.0625, fi c = 5.0148, fi p = 1.27973, L = 9 3 .) The legend describes the numbers of 
occurrence for each data point (N c , N p ). For better visibility data are grouped in bands. Note 
that x- and y-axis have different scales. P(N c ,N p ) can be described as a bivariate Gaussian in 
N c — (N c ) and N p — (N p ). Susceptibilities XTp p ari ^ Xt°p c can ^ e extracted from the half- widths 
of these distributions for p p =)p p (= const and p c = (p c ) = const (small black bars), respectively. 
Similarly, Xt jl p an d Xt^ can ' 3e obtained from the half-widths of the projections to the x- and y- 
axis (large black bars). We can also define an order parameter along the main axis of the ellipsoidal 
contours which will maximize fluctuations and result in the order parameter susceptibility x+ as 
described in the text. 
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FIG. 13: Intermediate incoherent structure factor of colloids and polymers plotted versus time (note 
the logarithmic scale) for rj p = 0.065,0.129,0.197,0.255,0.318 and the wave vectors (a) q = 0.93 
and (b) q = 6.1. For F£(q,t) only one curve (dashed) is shown as it hardly changes with polymer 
concentration. 
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FIG. 14: Plot of (1/6) dg a (t)/dt versus t (note the logarithmic scale) for colloids (upper part) and 
polymers (lower part) for the same choices of ry p as in Fig. [T3j 
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FIG. 16: Mean square displacement, Eq. (f4T|) . relating to interdiffusion (upper part) and its time 
derivative A(t), Eq. (fl2l) . (lower part). The data shown corresponds to the same values of e as in 
Fig. [ID 
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FIG. 17: Onsager coefficient A for interdiffusion plotted versus e = 1 — Typ/^p^rit (symbols with 
error bars) . Full circles show the prediction of the Darken equation (|45l) . 
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